Dynamical Preeze-out in 3-Fluid Hydrodynamics 
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Freeze-out procedure accepted in the model of 3-fluid dynamics (3FD) [J] is analyzed. This 
procedure is formulated in terms of drain terms in hydrodynamic equations. Dynamics of the 
freeze-out is illustrated by 1-dimensional simulations. It is demonstrated that the resulting freeze- 
out reveals a nontrivial dynamics depending on initial conditions in the expanding "fireball". The 
freeze-out front is not defined just "geometrically" on the condition of the freeze-out criterion met 
but rather is a subject the fiuid evolution. It competes with the fluid flow and not always reaches 
the place where the freeze-out criterion is met. Dynamics of the freeze-out in 3D simulations is 
analyzed. It is demonstrated that the late stage of central nuclear collisions at top SPS energies is 
of the form of three (two baryon-rich and one baryon-free) fireballs separated from each other. 
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I. INTRODUCTION 

Hydrodynamics is now a conventional approach to sim- 
ulations of heavy- ion collisions. Even review papers 
@j B S H, Hi do not comprise a complete list of numerous 
applications of this approach. The hydrodynamics is ap- 
plicable to description of hot and dense stage of nuclear 
matter, when the mean free path is well shorter than the 
size of the system. However, as expansion proceeds, the 
system becomes dilute, the mean free path becomes com- 
parable to the system size, and hence the hydrodynamic 
calculation should be stopped at some instant. All hy- 
drodynamic calculations are terminated by a freeze-out 
procedure, while these freeze-out prescriptions are some- 
what different in different models. 

In the present paper we would like to describe in more 
detail the freeze-outprocedure accepted in recently de- 
veloped 3FD model [l|, 0] • The 3FD model is designed for 
simulating heavy-ion collisions in the energy range from 
BNL Alternating Gradient Synchrotron (AGS) to CERN 
Super Proton Synchrotron (SPS). Unlike the conven- 
tional hydrodynamics, where local instantaneous stop- 
ping of projectile and target matter is assumed, a spe- 
cific feature of the dynamic 3-fluid description is a finite 
stopping power resulting in a counter-streaming regime 
of leading baryon-rich matter. The basic idea of a 3-fluid 
approximation to heavy-ion collisions [l|,0] is that at each 
space-time point a generally nonequilibrium distribution 
function of baryon-rich matter, can be represented as a 
sum of two distinct locally equilibrated contributions, ini- 
tially associated with constituent nucleons of the projec- 
tile (p) and target (t) nuclei. In addition, newly produced 
particles, populating the mid-rapidity region, arc associ- 
ated with a "flreball" (f) fluid. This model is a straight- 
forward extension of the 2-fluid model with radiation of 



direct pions 0, i, and (2-hl)-fluid model 
In particular, the 3FD model allows a certain formation 
time for the fireball-fluid production, during which the 
matter of the fluid propagates without interactions. 
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We have started our simulations P, 
pic hadronic equation of state (EoS) 
natural reference point for any other more elaborate EoS. 
The 3FD model turned out to be able to reasonably re- 
produce a large body of experimental data [l], [l^, [ij] in a 
wide energy range from AGS to SPS. This was done with 
the unique set of model parameters summarized in Ref. 
IL Problems were met in description of transverse flow 
ISj . The directed flow required a softer EoS at top AGS 
and SPS energies (in particular, this desired softening 
may signal occurrence of a transition into quark-gluon 
phase) . 

In particular, transverse-mass spectra of various 
hadrons were reproduced [11, Experimental data 

on transverse-mass spectra of kaons produced in central 
Au+Au or Pb-hPb [l^ collisions reveal peculiar de- 
pendence on the incident energy. The inverse-slope pa- 
rameter (so called effective temperature) of these spectra 
at mid rapidity increases with incident energy in the AGS 
energy domain and then saturates at the SPS energies. 
In Refs. [H, [l^ it was assumed that this saturation is 
associated with the deconfinement phase transition. This 
assumption was indirectly confirmed by the fact that mi- 
croscopic transport models, based on hadronic degrees 
of freedom, failed to reproduce the observed behavior of 
the kaon inverse slope po| . Hydrodynamic simulations 
of Ref. [U succeeded to describe this behavior provided 
the incident-energy dependence of the freeze-out temper- 
ature has a very similar shape to that of the correspond- 
ing kaon effective temperature. Thus, the puzzle of kaon 
effective temperatures was just translated into a puzzle 
of freeze-out temperatures. 

In Ref. [31 it was shown that dynamical description 
of freeze-out, accepted in the 3FD model, naturally ex- 
plains the incident energy behavior of inverse-slope pa- 
rameters of transverse-mass spectra observed in experi- 
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ment. This freeze-out dynamics, effectively resulting in a 
pattern similar to that of the dynamic liquid-gas transi- 
tion, differs from conventionally used freeze-out schemes. 
This is the prime reason why we would like to return 
to discussion of assumptions underlying this prescription 
and present one-dimensional simulations, clarifying con- 
sequences of this freeze-out. It is natural to start this 
discussion with a critical review of standard assumptions 
of the freeze-out and recent developments in this field. 

II. FREEZE-OUT: STILL DEBATED PROBLEM 

The hydrodynamic simulation is terminated by a 
freeze-out procedure. Though this method (as applied to 
high-energy physics) was first proposed almost 50 years 
ago by Milekhin [24I, this is a still debated problem. 
The method was intuitively clear and easily applicable. 
However, Cooper and Frye [l^l claimed that Milekhin's 
method violates the energy conservation. To remedy the 
situation, they proposed their own recipe, in which the 
observable spectrum of a-hadrons is calculated as follows 

e'^^ jja{p^K) Up,x) (1) 

where E is a 3-dimensional hypersurfacc on which a cer- 
tain criterion of the freeze-out is met. Here integration 
runs over this hypersurfacc, n'^ is normal vector to the 
element da of this hypersurfacc, and fa{p,x) is an equi- 
librium distribution function of a-hadrons 

^) = (1^ CXp{{p,uf^^f,a)/T}±l 

defined in terms of local thermodynamic and hydrody- 
namic quantities on this freeze-out hypersurfacc: chem- 
ical potential /ia(a:), temperature T{x) and 4- velocity 
u^{x). Here ga is degenerac y o f the a particle. 

The Cooper-Frye recipe [23| is now extensively used 
in hydrodynamic calculations, see, e.g., [l^, [H, [H, [13, 
m. Hi, [lo, im, m, m, m, lll, m. However, it is not 
free of problems neither. It gives negative contribution 
to the particle spectrum in some kinematic regions in 
which the normal vector to the freeze-out hypersurfacc 
is space-like, p^n'^ < 0. This negative contribution cor- 
responds to frozen out particles returning to the hydro 
phase. Cut off of this negative contribution again re- 
turns us to the violation of the energy conservation. To 
get rid of this negative spectrum, there was proposed a 
modification of the Cooper-Frye recipe based on a cut- 
Jiittner distribution [stI . IssI . ISa . |40| . In this distribution 
the part of the Jiittner distribution that gave the negative 
spectrum is simply cut off. To preserve the particle and 
energy conservation, the rest of Jiittner distribution is 
renormalized, effectively resulting in a new temperature 
and chemical potential (so called "freeze-out shock"). In 
fact, this cut-Jiittner recipe has no physical justification, 
except for practical utility. Moreover, the cut-Jiittner 



recipe is not supported by schematic kinetic treatment 
[4H of the transitional region from hydro regime to that 
of dilute gas. Recently there was proposed a new freeze- 
out recipe, a canceling-Jiittner distribution [4^ . which 
complies with results of schematic kinetic treatment [4l| . 
It should be stressed that this was precisely the schematic 
kinetic treatment. This region, where the transition from 
highly coUisional dynamics to the coUisionless one occurs, 
is highly difficult for the kinetic treatment and hardly al- 
lows any justified simplifications. 

All above considerations of the freeze-out process, in- 
cluding both the original Cooper-Frye prescription and 
its improvements, proceed from the following assump- 
tions: 

I. "Decoupling" of matter from hydrodynamic regime 
happens on a continuous hypersurfacc S. 

II. This hypersurfacc is determined on the requirement 
that a certain criterion of the freeze-out is met: 
e.g., temperature, energy density or baryon density 
reaches a certain value. 

HI. After this "decoupling" particles stream freely to 
detectors. 

In fact, transition from highly coUisional (hydro) regime 
to coUisionless one occurs in some finite 4-volume. As- 
sumption (I) is just an idealization — this 4-volume is 
shrunk to a hypersurfacc. Conservation conditions on 
such hypersurfacc are constructed in analogy with shock 
front in hydrodynamics and result in the Cooper-Frye 
formula ([T]). However, the requirement that this surface 
is continuous does not follow from anywhere. It is just 
an assumption. For instance, if we assume a discontin- 
uous hypersurfacc, i.e. that consisting of tiny (infinitely 
small in continuum limit) fragments with normal vectors 
coinciding with local 4- velocity, nj^ = ^ then we return 
to the original Milekhin's method of the freeze-out. 

The Milekhin's method assumes that a hydro system 
freezes out by emitting tiny fireballs of matter. Let P^qj 
is the total 4-momentum of the system. Then at the 
first step of the freeze-out a tiny droplet with the 4- 
momentum AP/* is emitted 

PL^PU + ^Pt, (3) 

where P^^j^ is the 4-momentum of the still hydro- 
evolving fiuid. In terms of the energy-momentum tensor 

T/'!' . the AP/" 4-momentum can be written out as follows 

I*; ' ' 

^Pt=l dVT^^^l da^nl, (4) 

where AK; is the volume of the fireball in the reference 
frame, where T^^^ is considered. The last equality in Eq. 

(HI represents AP/" in the covariant way, i.e. in terms 
a hypersurfacc element AE^ and the normal vector to 
this element n^, cf. Eq. ([T|). In particular. Milekhin's 
choice consists in rij^ = m^. From representation ^ it 



3 



may seem that relation between AP/' and T^^'^ depends 
on AEi. This would imply that a proper hypersurface 
element AE^ should be chosen to maintain relation 
In fact, the r.h.s. of Eq. (|4]) is independent of n't. The 
formal proof of that can be found, e.g., in Ref. [4^. It 
is possible to demonstrate it in a simpler way. Let us 
write down Tf^? in terms of contributions of individual 



particles [43 



(5) 



where p!^ (t) and x„ are the 4-momentum and the instant 
position of the nth particle, respectively. Integrating ex- 
pression ([5]) over volume AVi, accordingly to Eq. we 
arrive at 



(6) 



Here spurious dependence on AE; reveals itself in a seem- 
ing dependence of the r.h.s. of Eq. ([6]) on the syn- 
chronized time instant t, which really depends on the 
reference frame and hence on AE,;. Note that the Pj^ 
quantity is assumed to be conserved, therefore the time 
dependence of the r.h.s. of Eq. ([6]) is completely inap- 
propriate. 

Let us consider the r.h.s. of Eq. ^ in two refer- 
ence frames, i.e. on two hypersurface elements AE^ and 
AE^. The time synchronization depends on the reference 
frame. Therefore, in the sums over particles 



and 



(7) 



some pl^{t) and p'n{t') may occur, which are not sim- 
ply related by the Lorentz transformation but are com- 
pletely different because the corresponding particles at 
the t' instant have exercised additional interactions (or 
vice versa, have not exercised all those interactions) as 
compared to those completed to the t instant. And nev- 
ertheless two sums in Eq. ([7]) are equal, since in each 
two-particle or multi-particle point-like interaction the 4- 
momentum is conserved. The point-like character^ of 
the interaction is of prime importance here. If particles 
interact point-like, they change their momenta simulta- 
neously in any reference frame. Thus, the r.h.s. of Eq. 
(ini) is really independent of time t and hence of AE^. 

In view of Eqs. ^ and upon completion of the 
freeze-out process, we have 



(8) 



^ Action at distance in the relativistic case requires introduction 
of fields mediating this interaction. Then the field contribution 
should be also included in T^^ . For the sake of simplicity, we 
confine ourselves to the point-like interaction. 



where the hypersurface E consists of elements AE^. As 
we have seen, this 4-momentum conservation does not 
depend on the choice of hypersurface elements AE^. 



Milekhin's choice is 



and results in a discon- 



tinuous hypersurface. The Cooper-Frye choice proceeds 
from requirement of continuity of the E hypersurface. 
Difference between these two choices is illustrated in Fig. 
[T] The lower panel of Fig. [1] shows a schematic structure 
of Milekhin's hypersurface. In practical calculations the 
fragments of Milekhin's hypersurface are so tiny that the 
whole hypersurface looks like in upper panel of Fig. [U 
however, with normal vector to each tiny fragment coin- 
ciding with the 4-velocity. 




FIG. 1: (Color online) Freeze-out hypersurface for hydrody- 
namic evolution the ID step-like slab of nuclear matter (see 
subsect. nil C|l . Initial conditions for this slab are constructed 
on the assumption that they are formed by the shock-wave 
mechanism in head-on collisions of two ID slabs at -Ei^b = 10 
A GeV. The upper panel displays the Cooper-Frye choice for 
the hypersurface. The lower panel schematically illustrates 
Milekhin's choice for the hypersurface. Arrows indicate local 
4-velocities on this hypersurface. 

Therefore, Milekhin's method in fact conserves the en- 
ergy, but to see it one should consider it on a discon- 
tinuous hypersurface. The baryon number conservation 
can be demonstrated in a similar way. The fact that 
different single-particle distributions (i.e. Cooper-Frye, 
cut-Jiittner with renormalization, canceling-Jiittner, and 
even Milekhin's distributions) provide all the required 
conservation laws and at the same time produce differ- 
ent particle spectra, cf. Eq. ^ and Refs. [H, SI], only 
increases ambiguity of freeze-out consequences and indi- 
cate the barest necessity for further studies of the freeze- 
out. Such studies of the freeze-out in a finite 4- volume 
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(a finite space-like layer [45[) are now in progress. 

From the practical point of view, assumption (II) 
means that we should first run the hydro calculation 
without any freeze-out and only after that look for a hy- 
persurface, where the freeze-out criterion is met. This hy- 
persurface is determined as if the hydrodynamic system is 
not affected by the freeze-out. This procedure indeed re- 
sults in a continuous hypersurface, which could justify as- 
sumption (I). The method of "continuous emission" [i^ 
offers a more consistent way of performing freeze-out. 
This method considers a continuous emission of particles 
from a finite volume, governed by their mean free paths. 
In this approach the freeze-out process looks like an evap- 
oration (or fragmentation, on account of final-size grid) 
first from the system surface and then as a volume frag- 
mentation of the system residue. The particle emission 
from surface layer of the mcan-free-path width may hint 
at a discreetness of the hypersurface arising when this 
width is shrunk to zero. 

In particular, the continuous-emission mechanism im- 
plies that dynamics of this evaporation, characterized 
by its own rate, competes with the hydrodynamic ex- 
pansion. Therefore, freeze-out front may not reach the 
places, where, e.g., the energy density reaches the criti- 
cal value. This is already a consequence of real freeze-out 
dynamics. Unfortunately, this method is very difficult for 
the numerical implementation because probability for a 
particle to leave the hydrodynamic system depends not 
only on the past but also on the future evolution of this 
system, since particle emission occurs from time-evolving 
system. 

Assumption (III) is also an approximation to reality. 
This was the reason why cascading was applie d after the 
hydrodynamic freeze-out in Refs. [26l. [29l. Is^. Issj . This 
cascading allowed, in particular, to reproduce a two-slope 
form of transverse-mass spectra [1^ and a correct value 
of elliptic flow [23, HI] . Inclusion of some inelastic chan- 
nels in this cascading may be important for proper re- 
production of multiplicities, e.g., the K~ multiplicity 
However, even this cascading is not enough. A mean- 
field cascading is really needed. The reasons for this are 
as follows. First of all, the matter at the freeze-out in- 
stant is still dense enough, such that a part of energy is 
accumulated in collective mean fields. This mean-field 
energy should be released before calculating observables. 
In the presently discussed 3FD model [l| we do this at the 
freeze-out stage by recalculating thermodynamic quanti- 
ties in terms the hadronic gas EoS rather than a nontriv- 
ial EoS used in the hydro computation. 



III. FREEZE-OUT IN 3FD 

The freeze-out scheme adopted in the 3FD model 
is an attempt to modify and improve the standard 
freeze-out procedure in certain aspects rather than a 
final solution of the freeze-out problem. Let us start with 
criteria of the freeze-out. We formulate these criteria in 



terms of energy density, which is a universal quantity 
applicable both at very high energies (instead of tem- 
perature) and at low energies (instead of baryon density). 

(i) The freeze-out criterion we use is 

e < Efrz, (9) 

where 

£ = u^r^'u, (10) 

is the total energy density of all three fluids in the proper 
reference frame, where the composed matter is at rest. 
This total energy density is defined in terms of the total 
energy-momentum tensor 

T^"" = T^" + rr + 7/^" (11) 

being the sum over energy-momentum tensors of sepa- 
rate fluids, and the total collective 4- velocity of the mat- 
ter 

= u^Tf^yiuxT^^u,). (12) 

Note that definition is, in fact, an equation deter- 
mining M^. In general, this does not coincide with 4- 
velocities of separate fiuids. This definition of the collec- 
tive 4-velocity is in the spirit of the Landau-Lifshitz ap- 
proach to viscous relativistic hydrodynamics [47j . Only 
the formed (to the time instant of consideration) part of 
the f-fluid is taken into account in T/"' of Eq. pTjl . To 
the end of the freeze-out process all f-fluid turns out to 
be formed and hence frozen out. 

In the present simulations wc use the value Sf^z — 0.4 
GeV/fm'^ as the critical freeze-out energy density, with 
the exception of low incident energies i?iab, for which 
we use lower values: efrz{'2A GeV) = 0.3 GcV/fm'^ and 
EfrzilA GeV) = 0.2 GeV/fm3. 

(ii) In order to prevent freeze-out of initial cold 
nuclei, we apply the additional criterion 

Ufj^d^e < at the system surface, (13) 

i.e. at the boarder of matter with vacuum, which coin- 
cides with the freeze-out front. In the frame, where the 
freeze-out front is at rest, dt£ = 0, condition reduces 
to uVe < 0, which demands that collective velocity of 
the matter is directed outside the system. To meet this 
condition, wc in fact start the freeze-out procedure only 
at the expansion stage of the collision (see next subsect.). 

(iii) A very important feature of our freeze-out 
procedure is an anti-bubble prescription, preventing 
formation of bubbles of frozen-out matter inside the 
dense matter still hydrodynamically evolving. The 
matter is allowed to be frozen out (provided two above 
criteria are met) only if 

(a) either the matter is located near the boarder with 
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vacuum (this piece of matter gets locally frozen out) 
(b) or the maximal value of the total energy density in 
the system is less than eftz 

maxe < Sfrz (14) 

(the whole system gets instantly frozen out). 



where symbol (s) in the subscript indicate that the value 
is taken at the matter side of the freeze-out front. The 
latter is the effect of S +0 introduced in Eq. We 
have also omitted friction forces in the r.h.s. of Eq. (|16p . 
since they are really unimportant at the freeze-out stage. 
Integrating above equations over small Aa; around a; = 0, 
we arrive at 



Criterion (iiib) is convenient for numerical implemen- 
tations while docs not look quite physical. From physi- 
cal point of view, it would be preferable to change max e 
to the energy density averaged over the system, (e). In 
view of the discussion below, such a substitution will not 
change the qualitative pattern of the freeze-out, however, 
can somewhat affect quantitative results at AGS energies. 

Before the instant of the global freeze-out, cf. (iiib), 
above freeze-out criteria can be summarized in terms of 
dynamic equations 
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d^^T^ = (Friction)" 



(15) 



(16) 



where Jj^ is the baryon current of the a fluid, a = p, t 
or f (i.e. projectile, target or fireball), note that Jf* = 0. 
Here (Friction)'' stands for interaction terms between flu- 
ids, the explicit form of which is not important here. Og 
is a step function at the system surface, which takes into 
account criterion (iiia) . 1-dimensional simulations of the 
freeze-out (i)-(iiia) show that it results in discontinuity 
of e (and other quantities) at the system surface. This 
discontinuity is numerically smeared out only to the ex- 
tent of the finite grid step. Therefore, in analytic equa- 
tions the step function Qs can be represented by a sharp 
step function 



e, = e(e - s) 



with S 



-0. The function 



(17) 



(18) 



with e" being the value of e at matter side of the surface 
discontinuity, takes into account conditions (i) and (ii) 
of the freeze-out. 

Let freeze-out conditions (i) and (ii) be met, i.e. 
Bfrz = 1. To verify that Eqs. (IT5l)-((T8l) actually cor- 
respond to the above described scheme (i)-(iiia), let 
us consider these equations for the stationary situation 
in the reference frame, where the freeze-out front is at 
rest, i.e. dtQg = 0. For the sake of convenience, let 
us associate the freeze-out front with the x = plane. 
Then the only nonzero component of the 9^ 8s 4- vector 
is dx&s = (the hydrodynamic matter occupies the 

a; < semi-space), and Eqs. (fT5 l) - ([T6l) take the form 



'Jx'Ja 



ix>0) 

rpXV 

^a(x>0) 



rx TX 



rjiXU 



rpXV 



-TZLMx), 



Thus, terms oc i9^8s in the r.h.s. of Eqs. and m 
play role of sinks, which remove matter from hydrody- 
namic evolution, making the hydrodynamic quantities 
>^o!(x>o) and T„(a;>o) to be zero after the freeze-out front. 

This kind of freeze-out is similar to the model of "con- 
tinuous emission" proposed in Ref . (46j . There the parti- 
cle emission occurs from a surface layer of the mean-free- 
path width. In our case the physical pattern is similar, 
only the mean free path is shrunk to zero. 

The above discussion concerns only the first part of the 
freeze-out procedure, i.e. the application of freeze-out 
conditions. The second part consists in calculation of 
spectra of observable particles: 

(iv) Milckhin's method [l^, defined on a discontinuous 
hypersurface was used for calculation of observables, cf. 

sect, ini 

Thus, this freeze-out is similar to evaporation of par- 
ticles from a free surface of the system followed by ex- 
plosion of the fluid residue, if criterion ^ is met in the 
whole finite volume of this residue. 



A. Physical pattern of the freeze-out 

The physical pattern behind this freeze-out resembles 
the process of expansion of compressed and heated clas- 
sical fluid into vacuum. Physics of this process is studied 
both experimentally and theoretically [H, |4^, \E^, [5l| . 
Evaporation from free surface of normal (not super- 
heated) fluid is a very slow process. Accordingly, the 
freeze-out of matter of high density {e > £frz) is sup- 
pressed in our model. During expansion the fluid be- 
comes more rarefied, still remaining quite hot. Thus, the 
fluid becomes superheated at e < e^z- It occurs first at 
the periphery of the system, which is first affected by the 
decompression wave. Evaporation from free surface of 
superheated fluid is already a fast process. Accordingly, 
the freeze-out is allowed at e < etrz- 

Situations are possible, when freeze-out criterion ^ is 
met in the whole slab near the free surface rather than 
only at the surface. Such situations are illustrated in 
Subsects. nil CI and IIII Dl Here we have a choice either 
to instantaneously freeze out this whole near-surface slab 
or to wait until the freeze-out front will gradually traverse 
this slab (if ever). Making this choice, we rely on results 
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of experiments on evaporation from superheated fluids. 
It was shown (see, e.g., Ref. [5l|) that the evaporation 
front propagates with respect to fluid not faster than with 
the speed of sound. Precisely this choice is realized in our 
model by means of condition (iiia) or alternatively by 
dynamic equations (|15|) and ()16|1 . Thus, the freeze-out 
front may stay at essentially lower energy densities than 
Efrz because supersonic fluid expansion prevents it from 
reaching the region, where e = et^z. 

Physically it implies that a particle is evaporated 
("frozen-out") only if it escapes from the system with- 
out collisions. Thus, its mean free path (Amfp) should be 
larger than its path to the free surface (with due account 
of the future evolution of the fluid). Precisely this cri- 
terion is applied in the model of "continuous emission" 
(46| . In our simplified version of the "continuous emis- 
sion" , Amfp = in the fluid phase and Amfp oo in the 
gas phase. Therefore, a particle can escape only from the 
free surface, which cannot move inward the system faster 
than with the speed of sound [5l|. 

The only exception from this rule we do at the final 
stage of the freeze-out. As it was observed in experiments 
with classical fluids (see, e.g., Ref. [i^), a fluid trans- 
forms into gas by explosion, if it is strongly superheated 
all over its volume. Therefore, at the final stage of the 
freeze-out, when criterion ^ is met in the whole volume 
of the fluid residue, we assume that the whole residue 
becomes frozen out simultaneously (condition (iiib)). 

Of course, criterion ^ is not universal. In particular, 
it is not applicable to the cold nuclear matter, which 
has £ « 0.15 GcV/iia^ in its ground state. Additional 
condition (jl3p . preserving the cold nuclear matted from 
being frozen out, is directly connected with this fact. We 
hope that criterion ([9]) is good enough for a restricted 
domain of the phase diagram, where freeze-out of hot 
nuclear matter really occurs. 

In order to further clarify our freeze-out scheme it is 
useful to consider the way it was implemented in the 
numerical scheme. 



B. Numerical implementation of the freeze-out 

The numeric scheme of the 3FD code is based on the 
modified particle-in-cell method [s^ . [ssj , which is an ex- 
tension of the scheme first applied in Los- Alamos . In 
the particle-in-cell method each time step of the compu- 
tation consists of three stages: (I) Eulerian stage, (II) La- 
grangian stage, and (III) transformation from the frame 
of computation to the local rest frame of fluids in order 
to calculate thermodynamic quantities and flow veloci- 
ties (as applied to the freeze-out, it is described in point 
(c) of this subscct.). 

The transfer of energy-momentum due to pressure gra- 
dients, friction between fluids and production of the fire- 
ball fluid is computed on the spatially flxcd grid (so called 
Eulerian stage of the scheme). The convective transfer of 
the baryonic charge, energy and momentum is performed 



at the Lagrangian stage of the scheme. At this stage the 
matter is represented by an ensemble of Lagrangian par- 
ticles which accumulate all the energy, momentum and 
baryon charge of the system. At the time step (let it be 
"1"), when the freeze-out has not started yet, the total 
energy-momentum (P/^j) and baryon charge (-Btot) are 
presented by sums over these particles 

P^ot = ^AP^(ti)-^Ay,„(ti)T(':°)(ii), (19) 

Stot = ^AB,„(ti) = ^Al/„(ti) 4„)(ti), (20) 

ia ia 

where Ai-*^(ti) and /S.Bia,{ti) are respectively the 
energy-momentum and baryon charge of an ith parti- 
cle belonging to the a-fluid. These sums run over both 
formed and still unformed particles of the f-fluid. In fact, 
baryon charges of all particles belonging to a baryon-rich 
fluid are taken to be constant and equal, lS.Bip{t) = bp 
and ABit{t) = bt {bp ~ bt if nuclei are identical), while 
for the f-fluid ABif{t) = 0. Each Lagrangian particle is 
also characterized by a volume AVia{ti). Therefore, the 
above quantities of the particle are expressed in terms 
of respective densities T/^'^x and J?. as it is indicated 
in Eqs. and (PO)) . Simulation is performed in the 

frame of equal velocities of colliding nuclei. Hence, all 
the quantities in Eqs. (|19p and ([^0]) are related to this 
frame. Eqs. and (PO]) imply that all the matter 

participates in dynamical evolution before the freeze-out 
starts: 

Pdynih) = PL, i?dy„(il)=i3tot. (21) 

Here P^y^ and B^yn are respectively total energy- 
momentum and baryon charge participating in dynam- 
ical evolution. We avoid term "hydrodynamic evolution" 
because a part of the f-fluid may be still unformed. 

In the present scheme the Lagrangian particle has 
a proflle function of the form and size of the grid cell 
with uniform distribution of densities. Therefore, in 3D 
simulations a single Lagrangian particle contributes to 8 
cells on the grid, with which it overlaps. These spatially 
extended particles make the scheme smoother and hence 
more stable. In Refs. P, this numerical scheme is 
described in more detail. 

(a) To roughly meet condition (ii), the freeze-out 
procedure is started only at the expansion stage of the 
collision, i.e. after the time interval required for nuclei 
to traverse each other, provided they keep their initial 
velocities. In the cm. frame of two identical nuclei 
this time delay is At^z = Da/ (jc-mVcm), where Da/Tchi 
is the Lorcntz contracted diameter of the nucleus. 7cm 
and Vcm are the 7 factor and the initial velocity of the 
nucleus in the cm. frame, respectively. 

(b) The freeze-out criterion ^ is checked at the La- 
grangian stage of each time step (let it be nth step). In 
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the f-fluid only those particles are considered which have 
been formed to the time instant i„ . For each Lagrangian 
particle it is checked in all cells, which overlap with this 
considered particle, i.e. in 8 cells. 

If the freeze-out criterion is met in all these 8 cells and 
if at least one of these cells is "empty" (i.e. contains no 
centers of any Lagrangian particles), then this considered 
Lagrangian particle is counted as frozen out. This is re- 
alization of condition (iiia). This frozen-out Lagrangian 
particle is removed from further hydrodynamic evolution. 
By doing this we remove respective portions from sums 
^ and (HDD: 

PliyM ^ P^y^{tn-l)- E ^PtM, (22) 

ia frozen out at 

Bdyn{in) = Sdyn(iri-l) " ^ ABia{tn). (23) 

ia frozen out at tn 

Only parts Pl^y^-^{tn) and -Bdyn(^n) are kept in further 
dynamic evolution, unlike the conventional Cooper-Frye 
method. It is important to mention that had we kept this 
frozen-out particle in the dynamic evolution, its energy- 
momentum content would be changed at later time steps, 
and hence a part of its energy-momentum would be 
gained by other particles. Then we would face problems 
with the energy-momentum conservation, if we wanted 
to use the AP^^ quantity at the instant of its freeze-out 
for calculation of observable spectra. 

If the freeze-out criterion is met in all cells of the 
system, all Lagrangian particles are counted as frozen 
out at this time step, as it is required by condition 
(iiib). This is the end of the freeze-out process. 

(c) When the freeze-out process is over, i.e. P'^y,^ ~ 
and Bdyn — 0, we are left with ensemble of frozen-out 
Lagrangian particles which precisely obeys conservation 
laws: 

Ptot = APtajtia frozen-out), (24) 

ia 

Btot — ^ APiak^ia frozen-out)- (25) 

ia 

Here the summation runs over frozen-out (at different 
time instants ti^ frozen-out) particles, unlike sums (|19p and 
where this summation is associated with a fixed 
time instant. These frozen-out particles are precisely 
those droplets mentioned in Eq. ([3|). 

A frozen-out ia-particle is still characterized by five 
hydrodynamic quantities, J^^ and T!^°, and volume 
Al/ict, all these in the reference frame of computation. 
For the calculation of spectra we need distribution func- 
tion formulated in terms of thermodynamic quantities: 
temperature (T*"), baryon (/x^") and strange (/i!,") 
chemical potentials, and hydrodynamic 4- velocity (wfo,), 
cf. Eqs. H]) and In fact, this recalculation of the 
hydrodynamic quantities into thermodynamic ones is 
performed at each time step of the scheme based on the 



nongas EoS (involving some mean fields) accepted in 
the calculation. This EoS is not suitable for calculation 
of the spectrum of observable particles. First we should 
release the energy stored in mean fields. To do this, we 
calculate T-^-), ^;^^^'^\ ^.r^^"''^ and vf^^^^^^^ based on 
the hadronic gas EoS and proceeding from conservations 
of total energy-momentum, baryon and strange charges 
in the frozen-out particle. 

(d) At this stage we are still free to choose either Cooper- 
Fryc or Milckhin's scheme to calculate observables, as it 
was argued in Eqs. ©-([S]). We use Milekhin's method 
[2^ . defined on a discontinuous hypersurface consisting 
of tiny fragments coinciding with volumes of frozen-out 
Lagrangian particles AVia. The observable spectrum of 
hadrons is calculated as follows 

= E y^"""''P.<ai,.s^ /^"(Sas) ip) (26) 
ia 

where is the volume of the ia-particle in its rest 

frame, the sum runs over all frozen-out particles of all 
fluids, fia{ga.s){p) is tlic equilibrium distribution function 
defined already in terms of local gas thermodynamic 
(yia(gas)^ ^i«(gas) ^^^^ ^1"^^^^^^) and hydrodynamic 
("fa(gas)) quantities, cf. Eqs. (P) and 

In this prescription the baryon number and energy- 
momentum are precisely conserved by construction. It is 
worthwhile to mention that both the Cooper-Frye and 
Milckhin methods possess the same main problem: they 
both do not reject contributions of frozen-out hadrons re- 
turning into the hydrodynamic phase and do this to pre- 
cisely the same extent, see Appendix |X1 In particular, in 
our calculation this problem probably reveals itself in the 
failure of reproduction of the pion directed flow [l], [l^ . 
The advantage of Milckhin's method is just practical: 
with the exception of the pion directed flow it quite suc- 
ccssfuUy works in the 3FD model. The canceling-Jiittner 
recipe [421 overcomes the returning-hadrons problem of 
the above methods. It would be of interest to apply it in 
the 3FD model. 



C. One-dimensional simulations 

In order to clarify physics described by Eqs. (|15p - p^ 
let us consider ID simulations based on them. In Fig. [2] 
decays of step-like slabs of nuclear matter arc presented. 
These simulations have been performed till the time in- 
stant of the global freeze-out, i.e. when criterion (iiib) 
is met. The same EoS as that used in 3D simulations 
[Hi [13, [3 is accepted in the present calculations. First of 
all, we see that the freeze-out front is really step-like. It is 
smeared only over two cells (independently of their size) 
due to numeric scheme. Note also that for this step-like 
initial geometry the supersonic flow of matter (vx > Cg, 
where Cs is the speed of sound) always occurs beyond the 
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initial boundary of the slab, while the flow within the 
initial boundary is always subsonic {v^ < Cg). This is 
important for understanding results displayed in Fig. [21 




FIG. 2: (Color online) Evolution of the energy density during 
the decay of the ID step-like slab of nuclear matter. Solid 
lines display calculations without freeze-out, while dashed 
lines, with freeze-out. Initial conditions for these slabs are 
constructed on the assumption that they are formed by the 
shock-wave mechanism in head-on collisions of two ID slabs 
a-t -Blab ~ W A GeV for upper panel and -Ei^b — 2 A GeV 
for lower panel. Subpanels show zoomed regions of the freeze- 
out front at the fixed time instant t = 3 fm/c in terms of the 
energy density in separate cells (displayed by dots). 

There are three important velocities in this problem: 
the hydrodynamic velocity of the matter (vx) at the po- 
sition, where the freeze-out front occurs, the speed of 
sound Cs , and the velocity of transfer of the constant 
value of e = 0.4 GeV/fm^: 



e{vet + const, t) = 0.4 GcV/fm 



(27) 



In fact, Eq. (|?7)) is the equation for the hydrodynamic 
characteristic curve related to the 0.4 GeV/fm-^ value of 
the energy density. The freeze-out front, as defined by 
Eq. (fT5)) and cannot propagate in the fluid medium 
faster than with the local speed of sound, like any per- 
turbation in the hydrodynamics. 



Different dynamic patterns of the freeze-out fronts dis- 
played in two panels of Fig. [2] are associated with dif- 
ferent relations between above three velocities. In the 
upper case of eo — 3 GeV/fm'^ we have > 0. Hence, 
the point, where the freeze-out criterion (|18p starts to 
be met, is transferred farther and farther from the initial 
system boundary, hence the system expands. In this re- 
gion the flow is supersonic: > c^. At the same time, 
the matter velocity with respect to the characteristic ve- 
locity is^ {vx — Ve)/{1 — VxVe) < Cs, i.e. it is less than 
the speed of sound. Therefore, the fluid flow does not 
carry the freeze-out front away from the point of e = 0.4 
GeV/fm'^, where the freeze-out criterion ([18]) is met. The 
frcezc-out front stays at this position, see the zoomed 
subpancl in the top panel of Fig. [51 Since the matter 
within this freeze-out front is distributed over the range 
of < e < 0.4 GcV/fm'^; an average value for the actual 
energy density of the frozen-out matter (let us denote it 
as Eout) is approximately half of the freeze-out jump, i.e. 
£out ~ efrz/2 = 0.2 GcV/fm-^ in this case, see Fig. [31 
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FIG. 3: (Color online) Instant values of actual energy density 
of the frozen-out matter for decays of nuclear-matter slabs 
with various initial conditions (labeled by initial energy den- 
sities £o)- Initial conditions are constructed on the assump- 
tion of the shock- wave mechanism in head-on collisions of two 
slabs at Bjab = 2, 4, 6, 10, 20 and 40 A GeV (from bottom 
left to top right). The time evolution is displayed till the time 
instant of the global freeze-out (iiib). 

In the bottom panel of Fig. [21 i.e. Eq — 1 GeV/fm'^, we 
have We < 0. Here, the point, where freeze-out criterion 
p8)) starts to be met, cannot be reached by the freeze- 
out front. Indeed, the freeze-out front is first formed 
beyond the initial boundary of the slab, i.e. in the re- 



Note that velocities are added relativistically. 
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gion of the supersonic flow of the matter. Since it cannot 
be faster than the sound, it cannot overcome the "su- 
personic barrier" in the back direction in order to reach 
the e = 0.4 GeV/fm'^ point. It stays at the position of 
the "supersonic barrier" and does not move. Therefore, 
the decaying and frcczing-out system does not expand. 
Moreover, the front stays at the position of low energy 
density, e ~ 0.2 GeV/fm'^, see the zoomed subpanel in 
the bottom panel of Fig. [2l Therefore, the actual en- 
ergy density of the frozen-out matter turns out to be 
lower than efiz/2 (as in above case), i.e. £out ~ (0-2 
GeV/fm3)/2 0.1 GeV/fm^, see Fig. d 

In Fig. [3] we see that Eout remains practically constant 
during the major period of the freeze-out. The steep fall 
of Sout just before the global freeze-out occurs because 
the velocity on the characteristic curve (P7|) changes its 
sign. Then the point, where e — 0.4 GcV/fm'^ is achieved, 
starts rapidly move inwards the system, and hence the 
freeze-out front remains at lower energy density. 

Thus, we see that the freeze-out is not inseparably 
linked with the freeze-out energy density efrz but can oc- 
cur at lower energy densities due to dynamical reasons. 
At low initial energy densities the freeze-out front stays at 
lower energy densities than eftz and hence the actual en- 
ergy density of the frozen-out matter Eout turns out to be 
lower than £frz/2, see Fig. [3] With the initial energy den- 
sity rising, the £out quantity grows and gradually reaches 
the efrz/2 value and than approximately saturates. 




x[fm] 

FIG. 4: (Color online) Characteristic curves for hydrody- 
namic evolution the ID step-like slab with initial conditions 
corresponding to the upper panel of Fig. [2] Solid lines dis- 
play characteristic curves, corresponding to e = 0.4 and 0.7 
GeV/fm'^, calculated with local freeze-out (iiia) but without 
global freeze-out (iiib). Dashed lines correspond to the cal- 
culation without freeze-out. 

A standard procedure of performing freeze-out, which 
is applied in the major part of hydrodynamic calcula- 
tions now, proceeds in the following way. The hydro 
calculation runs absolutely unrestricted. The freeze-out 
hypersurface is determined by analyzing the resulting 4- 



dimensional field of hydrodynamic quantities on the con- 
dition of the freeze-out criterion met. In our case this 
would be the characteristic curve of e = 0.4 GeV/fm'^ 
calculated without freeze-out, displayed in Fig. [H 

In the 3FD model, the freeze-out criterion is checked 
continuously during the simulation. If some parts of the 
hydro system meet all criteria ((i), (ii), and either (iiia) 
or (iiib)), they decouple from the hydro calculation. The 
frozen-out matter escapes from the system, removing all 
the energy and momentum accumulated in it, cf. Eq. 
([3]). Therefore, it produces no recoil to the rest of still 
hydrodynamic system. The removal of the matter affects 
the system evolution. This influence is illustrated in Fig. 
m The £ = 0.4 GeV/fm'^ characteristic curves calculated 
with and without freeze-out turn out to be different. At 
the same time the e = 0.7 GeV/fm^ characteristic curves, 
which lie quite deep inside the system, remain unaffected 
by the freeze-out. The freeze-out hypersurface (i.e. curve 
in l-fl dimensions) for this case is presented in Fig. [T] 
It differs from the corresponding characteristic curve be- 
cause of the global freeze-out which occurs at time instant 
t ~ 9 fm/c. 

It is of interest to compare particle spectra predicted 
by different models for freeze-out. In Fig. [S] rapidity 
distributions of various hadrons are demonstrated for 
three different prescriptions of the freeze-out: the above- 
described model of dynamical freeze-out (Dynamical FO) 
[sec Eqs. Cooper-Frye freeze-out on character- 

istic curve (P7|) [Cooper-Frye FO], and Milekhin freeze- 
out on characteristic curve (p7)) [Milekhin FO]. These 
spectra were calculated within ID hydrodynamic expan- 
sion of the slab of initial width of 4 fm and initial en- 
ergy density £o = 3 GeV/fm'^. Initial conditions for this 
slab are constructed on the assumption that it is formed 
formed by the shock-wave mechanism in head-on colli- 
sion of two ID slabs at E'lab — 10 A GeV. Only nucleons, 
deltas and pions are included in the EoS used in this cal- 
culation. Displayed rapidity distributions are normalized 
to unit aria transverse to the direction (x) of longitudinal 
expansion. Freeze-out surfaces for these three models of 
freeze-out are displayed in Fig. [1] for Dynamical FO and 
Fig. m (dashed curve for £ = 0.4 GeV/fm'^) for Cooper- 
Frye FO and Milekhin FO. 

It is difficult to directly compare the Dynamical FO 
and Cooper-Frye FO, since the Dynamical FO surface 
is inappropriate for the Cooper-Frye FO and vise versa. 
This is because the Dynamical FO requires the hydrody- 
namics to be modified by removing the frozen-out matter, 
while the Cooper-Frye FO needs the hydro calculation to 
be run absolutely unrestricted. This is the reason why 
we consider the Milekhin FO, which takes place precisely 
on the same characteristic curve p7)) as the Cooper-Frye 
FO. In particular, this is the reason why we compare only 
two models in Fig. [S] 

As seen from Fig. [SJ baryon rapidity distributions are 
quite close to each other within different freeze-out mod- 
els. The Dynamical-FO and Cooper-Frye-FO are even 
strikingly close. Apparently, this occurs because they 
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FIG. 6: (Color online) Partial contributions to rapidity spec- 
tra of baryons (upper panel) and thermal pions (lower panel) 
from space-like and time-like parts of the freeze-out sur- 
face. Spectra are calculated within two models of freeze-out: 
Cooper-Frye freeze-out on characteristic curve (|27p [Cooper- 
Frye FO], and Milekhin freeze-out on the same characteristic 
curve [Milekhin FO] . 
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FIG. 5: (Color online) Rapidity distributions of baryons (up- 
per panel) , thermal pions (middle panel) and A-isobars (lower 
panel) calculated within three models of freeze-out: the model 
of dynamical freeze-out [Dynamical FO], [see Eqs. (P)|- (|26p ]. 
Cooper-Frye freeze-out on characteristic curve (|27p [Cooper- 
Frye FO] , and Milekhin freeze-out on characteristic curve (|27p 
[Milekhin FO]. 



are confined by the baryon number conservation. This 
is not the case for thermal pions (i.e. without contribu- 
tion of A decays). Difference between their distributions 
within different models is quite spectacular. 

As it is seen from Fig. [HI this difference mainly 
originates from contributions from the time-like"^ parts 



of the freeze-out surface, i.e. from those parts, where 
the the Cooper-Frye normal vector n^p is space-like: 
ncF • ncF < 0. Difference between Cooper-Frye and 
Milekhin recipes is most strongly pronounced in this case. 
Note that the Cooper-Frye FO even reveals its generic 
deficiency — its " time- like" spectrum becomes negative at 
midrapidity. Physically it means that the midrapidity 
region is most abundantly populated by particles which 
have to be returned to the hydro phase. At the same 
time, contributions from the space-like parts of the sur- 
face, where ucf ■ "-C-F > 0, are quite similar within 
Cooper-Frye FO and Milekhin FO. Small diference be- 
tween the Cooper-Frye normal time-like vector n^p and 
and the hydrodynamic 4- velocity results is tiny difference 
in "time-like" spectra. 

Thus, spectra of newly produced particles are quite dif- 
ferent in different models of freeze-out. At the same time, 
the numerical baryon-number and energy conservation is 



^ We prefer to specify hypersurfaces by the character of space-time 
intervals within them, similarly to Refs. l37l |40|| and contrary to 



Refs. [H,!!! 
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better than 1% in all considered models. We still face 
the problem: which of these models is physically true. 
Unfortunately, above-mentioned experiments on expan- 
sion of compressed and heated classical fluid into vacuum 
[isl. [49I. [sol. l5ll| cannot answer this question. As we have 
seen, particle number conservation (baryon number, in 
Fig. O makes these models fairly close to each other. 
We would like only to mention that the dynamical FO 
model naturally explains [31 the incident energy behav- 
ior of inverse-slope parameters of transverse-mass spectra 
observed in experiment. However, this is not a physical 
justification of this model. The question still remains 
open. 
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FIG. 7: (Color online) The same as in Fig. [2] (top panel) but 
for different values of Ax/Ai. Only results with freeze-out 
are displayed. 
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Results of the above reported ID simulations are sta- 
ble with respect to the numeric procedure [l^l- They 
are quite insensitive to the size of the cell. The freeze- 
out front is smeared only over two cells independently of 
their size and hence indeed is step-like in continuum limit. 
Dependence on the most important numeric parameter — 
the ratio of the space-grid step to the time step Ax/ At — 
is displayed in Figs. [7] and [51 To reduce numerical diffu- 
sion, this ratio should be taken optimal. As it was found 
in 1-dimensional simulations of exactly solvable problems 
[5^, the optimal range of this ratio is 2.5 < Ax/ At < 5 
with the preferable Ax/ At ~ 3.5, minimizing the numeri- 
cal diffusion. As seen, within this range 2.5 < Ax/ At < 5 
the results are really stable. 



D. 3D simulations 



Condition (i) (or Eq. ([T5)) ) ensures only that the ac- 
tual freeze-out energy density, at which the freeze-out 
actually occurs, is less than Shz- Therefore, e^z can be 
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FIG. 8: (Color online) The same as in Fig. [3] (bottom panel) 
but for different values of Ax/ At. 
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FIG. 9: (Color online) Actual average freeze-out energy den- 
sity in central {b — 0) Pb-|-Pb collisions as a function of in- 
variant incident energy. 



called a "trigger" value of the freeze-out energy density. 
As it has been demonstrated in the previous subsection, 
a natural value of this actual freeze-out energy density 
is Eout ~ £"72, i.e. at the middle of the fall from to 
zero. To find out the actual value of eout, we have to 
analyze results of a particular simulation. In our previ- 
ous paper [l|] we have performed only a rough analysis 
of this kind. This is why in the main text of Ref. [l| 
we mentioned the value of approximately 0.2 GcV/fm'^ 
for Sout and in the appendix explained how the freeze- 
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FIG. 10: (Color online) Time evolution (in the xz plane of the center-of-mass frame) of the total energy density (cf. Eq. 
(IIOI) ) in the central Pb+Pb collision at -Eiab = 158^4 GeV. The light-colored outer hallo corresponds to the simulation without 
freeze-out and thus indicates the matter which has already been frozen out to the time instant t. Arrows indicate velocities of 
baryon-rich fluids: open arrows for the projectile-like fluid and filled arrows for the target-like fluid. 
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out actually proceeded^. Results of more comprehensive 
analysis for central (6 = 0) Pb+Pb collisions are pre- 
sented in Fig. [9l which shows the Sout value averaged 
over space-time evolution of the collision: (eout)- As 
seen, (Sout) reveals saturation at the SPS energies. This 
happens in spite of the fact that our freeze-out condi- 
tion involves only a single constant parameter et^z = 0.4 
GeV/fm'^, with the exception of low incident energies, for 
which we use lower values: efrz(2A GeV) = 0.3 GeV/fm^ 
and eft.z(lA GeV) = 0.2 GeV/fm^. 

The "step-like" behavior of (£out) is a consequence of 
the freeze-out dynamics which has already been illus- 
trated in Fig. [3l At low (AGS) incident energies, the 
energy density achieved at the border with vacuum, e*, 
is lower than Eftz. The surface freeze-out stays at this 
lower energy density up to the global freeze-out because 
the freeze-out front cannot overcome the supersonic bar- 
rier in the expanding matter, cf. the lower panel of 
Fig. [51 At these low energies, the value (eout) turns 
out to be low sensitive to the freeze-out parameter Sfrz- 
Only the global freeze-out (iiib) of the system remnant, 
which also contributes to (cout); produces weak sensitiv- 
ity to Eftz. The values efrz(2A GeV) = 0.3 GeV/fm^ and 
efiz(l^ GeV) = 0.2 GeV/fm^ were precisely chosen in 
order to reduce contribution of the global freeze-out to 

(Sout)- 

With the incident energy rise the energy density 
achieved at the border with vacuum gradually reaches 
the value of e^z and then even overshoot it. If the over- 
shoot happens, the system first expands without freeze- 
out. The freeze-out starts only when e'* drops to the 
value of etrz- Then the surface freeze-out occurs really at 
the value e** w Sfrz and thus the actual freeze-out energy 
density saturates at the value (eout) ~ efrz/2. 

In Fig. [TO] the time evolution of the total energy den- 
sity (cf. Eq. (|10p ) in the central Pb-|-Pb collision at 
^'lab = 158A GeV is displayed. The hght-colored outer 
hallo corresponds to the simulation without freeze-out 
and thus indicates the matter which has already been 
frozen out to the time instant t. First of all we see that 
already in the beginning of expansion stage (t = 2 fm/c) 
the baryon-rich fluids are mutually stopped and unified 
to a good extent, since their hydrodynamic velocities al- 
most coincide: arrows, originating from the same point, 
are almost equal if not merged. This is not the case 
for the baryon-free fluid, since its formation lasts till ap- 
proximately t = 4 fm/c (see Fig. 18 in Ref. [ij). At the 
late stage of the expansion {t > 10 fm/c) the baryon-rich 
and baryon-free fluids become even spatially separated: 
The middle region of the system, containing no arrows, 
is solely populated by the baryon-free fluid, and hence 
the baryon-rich matter falls into two disconnected pieces. 



Thus, at the late stage of the evolution the system effec- 
tively consists of three "fireballs" (two baryon-rich and 
one baryon-free). This is in contrast to the assumption 
of the statistical model ([sl, [H, [13, H^), where a single 
uniform "fireball" is considered. The baryon-free "fire- 
ball" becomes frozen out first: the displayed time instant 
t = 11 fm/c is almost the last, when this "fireball" still 
hydrodynamically evolves. Evolution of the two baryon- 
rich "fireballs" till the complete freeze-out is rather long, 
till t « 20 fm/c. This spatial separation of "fireballs" 
happens only at high incident energies i^iab > 40^ GeV. 

The freeze-out in the longitudinal direction proceeds 
accordingly to the ID pattern of Fig. [2| (upper panel). 
In the transverse direction the freeze-out front moves in- 
wards the system. This a combined effect of fast lon- 
gitudinal expansion and comparatively slow transverse 
motion of the system. This effect actually results in the 
"two-fireballs" structure at the latest stage. 
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FIG. 11: (Color online) Time evolution (in the center-of-mass 
frame) of the non-frozen-out baryon number (normalized to 
the total baryon number of the system Ap + At) in central 
Pb-fPb collisions at various incident energies. 

The occurrence of two spatially separated "fireballs" 
at the latest stage of the collision is the reason why the 
global freeze-out, cf. (iiib), does not occur in this sys- 
tem. At low (AGS) incident energies, the surface freeze- 
out, actually occurring at lower than efrz densities, takes 
place only as long as the energy density in the center of 
expanding system exceeds eftz- When this center den- 
sity drops below Eftz, the rest of the system instantly 
gets frozen out, as it is illustrated in Fig. [11] at the 
example of time evolution the baryon number still in- 
volved in the hydrodynamic phase. The abrupt fall of a 
curve corresponds to the global freeze-out. In fact, only 
this global freeze-out stage depends on the freeze-out pa- 
rameter efrz at low incident energies. With the incident 
energy rise this instantly frozen-out remnant becomes 
smaller and smaller until it completely disappears above 
40^1 GeV energy. This happens because outer freeze-out 
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fronts stay already at the e^z density, therefore the maxi- 
mum volume density can be only higher than Sfj-z . At the 
same time the inner freeze-out fronts overtake these two 
baryon-rich "fireballs" from behind and result in gradual 
surface freeze-out of them without any remnant, where 
criterion (iiib) can work. 



IV. CONCLUSIONS 

The described method of freeze-out can be called dy- 
namical, since the freeze-out process here is integrated 
into fluid dynamics through hydrodynamic equations 
(fT5)) - p^ . The freeze-out front is not deflned just "ge- 
ometrically" on the condition of the freeze-out criterion 
met but rather is a subject the fluid evolution. It com- 
petes with the fluid flow and not always reaches the place 
where the freeze-out criterion is met. 

This kind of freeze-out is similar to the model of "con- 
tinuous emission" proposed in Ref. [1^ . There the parti- 
cle emission occurs from a surface layer of the mean-free- 
path width. In our case the physical pattern is similar, 
only the mean free path is shrunk to zero. In particular, 
the fact that the freeze-out sometimes occurs at lower 
energy densities than that prescribed by the freeze-out 
criterion can be associated with the dependence on the 
future evolution of the fluid system in the model of "con- 
tinuous emission" . In that model the particle emission 
occurs from interior of the time-evolving system. There- 
fore, if the emitted particle is able to leave the system or 
not depends on the expansion rate of the system. Similar 
situations happen in our model. Sometimes the freeze- 
out criterion is met too deep inside the system such that 
rapid expansion of the system prevents the freeze-out 
front from reaching this place. 

We argue that the original Milekhin's method [l^] of 
calculation of observable particle spectra is energy con- 
serving, as well as the Cooper-Frye recipe [1^, provided 
it is considered on a discontinuous hypersurface, i.e. that 
consisting of tiny (infinitely small in continuum limit) 
fragments with normal vectors coinciding with the lo- 
cal hydrodynamic 4-vclocity. Moreover, we argue that 
Milekhin's approach on a discontinuous hypersurface is 
natural for the Lagrangian formulation of hydrodynam- 
ics. There are no principal objections against using such 
a fragmented hypersurface instead of the continuous one 
like in the Cooper-Frye method. Discreetness of parti- 
cle emission may hint at discontinuous character of the 
freeze-out hypersurface. 

Systematic studies based on hydrodynamic models in- 
dicated that there are two distinct freeze-out points (see, 
e.g., [5^, [g^): chemical and thermal ones. Contrary 
to these studies we do not need two distinct freeze-out 
points. Multiplicities and spectra are simultaneously and 
quite satisfactorily reproduced with the single freeze-out 
described above [1[. Of course, we could somewhat reflne 
reproduction of experimental data by introducing two 
distinct freeze-out points. However, introduction of an 



additional fitting parameter would be too high price for 
such a slight improvement. 
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APPENDIX A: PROBLEM OF FROZEN-OUT 
HADRONS RETURNING INTO 
HYDRODYNAMIC PHASE 

When frozen-out matter coexists with still hydrody- 
namically evolving fluid, both the conventional Cooper- 
Frye recipe [23j and the model presented here suffer from 
a problem of frozen-out particles returning to the hydro 
phase. For the Cooper-Frye method this problem was 
discussed in a number of papers [13, [H, [H, H^] ■ Here 
we would like to demonstrated that precisely the same 
(and precisely to the same extent) problem exists for the 
freeze-out model discussed in this paper. 

Let the fluid and the frozen-out matter be separated by 
a space surface with an external normal 3-vector n. Let 
this surface move with velocity . Then we can decom- 
pose the number of particles frozen out on the element 
of hypersurface da, associated with the above separation 
surface, as follows (cf. Eq. ([1])) 

dN = dN,sc.+dN,et., (Al) 

diVcsc. ^dan^ [ f{p,x) e([v - v,] • n), (A2) 

J Po 

dNrct. =dan^ [ fip,x) e([v, - v] • n), (A3) 

J Po 

where v = p/po is the 3-velocity of a frozen-out particle. 
Here dN^sc. can be identified with a number of particles 
escaping from the system, since they move outward the 
system faster than the free surface overtakes them. While 
dNj-et. is a number of particles returning to the fluid. In 
the covariant forrn, 8-functions can be written as (see, 
e.g., Refs. [33,[3i,[3i,[43) 

e([v-v,]-n) = e(p,,n^^J, (A4) 
e([v,-v].n) = e(-p,.n^^J, (A5) 

where n^p^ is the normal 4-vector to the Cooper-Frye 
continuous hypersurface. This is always so, indepen- 
dently of the choice of 71^ we use for calculation of 
spectra. Therefore, the fraction of returning particles, 
dNrct./dNcsc. is independent of whether we employ the 
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Cooper-Frye recipe with = nJ^p^ or the Milckhin's method with = u^. 
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